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Abstract 

The exclusive queueing process (EQP) incorporates tlie exclusion principle 
into classic queueing models. It can be interpreted as an exclusion process of 
variable system length. Here we extend previous studies of its phase diagram by 
identifying subphases which can be distinguished by the number of plateaus in the 
density profiles. Furthermore the influence of different update procedures (parallel, 
backward-ordered, continuous time) is determined. 

1 Introduction 

Queueing theory is one of the most important topics in the field of operations research 
EIHIO. It has a broad spectrum of applications ranging from telecommunications to 
traffic engineering and supply chains. One of the simplest queueing processes is the 
so-called M/M/1 model, where customers enter the system with probability a and leave 
the system with probability /? at one server. The current state of the M/M/1 queueing 
process is completely specified by the number of customers. The system converges to 
a stationary state with a finite number of customers when a < p whereas the number 
of waiting customers diverges for a > /3. 

A feature which seems to be important for pedestrian queues and other traffic ap- 
plications is the excluded-volume effect: pedestrians can proceed only when there is 
enough space in front of them ||4l. This is e.g. seen in queues at the check-in at airports 
where passengers have to move the luggage when moving forward. However, standard 
queueing models like the M/M/1 model neglect the excluded-volume effect, and do 
not have a spatial structure. Then the length L of the system is given by the number 
of waiting customers (if customers have unit length) and the density is constant in 
space (p - N/L - 1). 

The "exclusive queueing process" (EQP) was introduced in |l5]0 to investigate how 
the excluded volume affects queues. It is obtained by modifying the input procedure 
of the one-dimensional totally asymmetric simple exclusion process (TASEP) with the 
ordinary "open boundary condition." Customers are injected always at the end of the 
queue, called left end in the following, i.e. behind the last customer waiting (Fig. [TJ. 
This is in contrast to the usual open TASEP where the input always happens at the 
same site, irrespective of the occupation of the other sites. The output is not changed 
compared to the usual open TASEP: in both models customers are extracted at the right 
end (the server) which is fixed. 

The EQP is not the only variant of the TASEP on a dynamic lattice. An ear- 
Uer example is the dynamically extending exclusion process (DEEP) introduced in 
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Figure 1: Exclusive queueing process. 



Q [8l 121 (To) as a model for fungal growth. In contrast to the EQP, the DEEP has no 
mechanism for reducing the system length and therefore the length of the system is 
always diverging. Other possible biological applications are length-regulation of mi- 
crotubules I.12J and bacterial flagellar growth [13J. 

The state space of the EQP is the set of configurations of the customers. It is more 
precisely given as 

5 ={T= lTn---ri|LeZ>o,T,=0,l), (1) 

where r = 0, 1 for L = 0, 1, respectively. The state t - It^-i ■ ■ ■ ti e S \{0) corresponds 
to the customer configuration where each site j is occupied or empty according to 
Tj = 1 or 0, and L defines the length of the system. The symbol corresponds to the 
state "no customer in the system." We denote the number of customers (#(Ty = 1)) by 
for a given state t = It/^-i ■ ■ ■ ti. Each customer enters the system, hops and leaves 
the system as 

^ 1 with probability a, 

1 ■ ■ ■ 1 1 ■ ■ ■ with probability a, 

■ ■ ■ 10 ■ ■ ■ ■ ■ ■ 01 ■ ■ ■ with probabihty p, ^ ^ 

■ ■ ■ 1 — > ■ ■ ■ with probability /3. 

The local update rules Q are not sufficient to specify the dynamics fully. In addition 
the sequence in which the rules are applied to the sites or particles needs to be defined. 
Here we consider two discrete-time updates (parallel and backward sequential updates) 
and a continuous-time dynamics where the parameters a, /3 and p are transition rates 
(not probabilities). 

In |5( il4il . exact stationary states for the continuous-time case and the parallel- 
update case were constructed as matrix product states based on the known forms for 
the corresponding TASEPs with a fixed length ifTSlfTBl . The phase boundary between 
the convergent and divergent phases was found to be modified compared to the classical 
M/M/1 queue. In particular, for the convergence to the stationary states, the injection 
rate (or probability) a cannot be bigger than the maximal current of the TASEP, i.e. 
"the queue itself is a bottleneck" as well as the server. 

In the phase diagram was analyzed in more detail. The convergent and di- 
vergent phases are both further subdivided in two subphases analogous to the maximal 
current and high-density phases of the TASEP. Furthermore time dependent properties 
were investigated. However the asymptotic form of the velocity for the growth of L 
in the divergent phase was left as an open problem. In this article we will clarify this 
point, giving density profiles with help of Monte Carlo simulations. (In the case where 
the customer hopping is deterministic (p = 1), an exact "dynamical state" in matrix 
product form exists which enables us to rigorously derive the behavior of {L,} and (A^,) 
CEJ.) 
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This article is organized as follows. In Section|2] we define the EQPs with various 
updates in more detail. The phase diagrams are derived and their relation with clas- 
sical queueing processes is discussed. In Section [3], based on simulation results we 
characterize further subphases of the divergent phases according to the shapes of the 
density profiles. In Section |4] we investigate the EQPs on the critical line separating 
the divergent and convergent phases. Finally, we give a summary and conclusions of 
this article in Section |5] 



2 Update rules 



The TASEP is a prototypical model of stochastic interacting particle systems. It has 
been studied intensively in the last decades both from the view of nonequilibrium sta- 
tistical physics 1 19, 20, 21 , 4| and mathematics (see e.g. I.22J ). Similar to the TASEP 
with the ordinary open boundary condition ll23l one can study the EQP with differ- 
ent update schemes. In the next two subsections we define two discrete-time EQPs 
and determine their phase diagrams. These are divided into four phases according to 
their asymptotic lengths (convergent vs. divergent) and the parameter-dependence of 
the outflows. Then we consider special cases and limits including the continuous-time 
EQP. 
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Figure 2: The phase diagrams of various queueing processes: the parallel EQP (top- 
left), the continuous-time EQP (top-middle), the backward EQP (top-right), the rule 
184 CA with stochastic boundaries (bottom-left), the continuous-time M/M/1 (bottom- 
middle) and the discrete M/M/1 (bottom -right). The relations among the phase dia- 



grams can be seen by ( 15 i. For example the bottom-left is obtained by setting p - \ in 
the top-left. 



2.1 Parallel update 

hi IS [mini [H, the EQP with the parallel update rule, which we call parallel EQP 
shortly, was investigated. In the parallel dynamics all sites are updated simultaneously. 
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e.g. 



101011 ^ 1011010 with probability ax /:> X (1 -p) Xy8, 
111011 ^ 110111 withprobability(l-Q')xpx(l-yS), (3) 
— > with probability I - a. 

The current-density relation for the TASEP with the parallel update and the periodic 
boundary condition is given by ||24l |25) 

which is also true for the ordinary open boundary condition lfT6l . The "critical line" 
that separates the parameter space into divergent and convergent phases is given by 
a - /||(pii), where 

•^"(Pii) = \ i_VI^ where A = W ^ m ^ r ^ ^^ 

with = 1 - -y/1 - p. When a < Juipn), the system converges to a stationary state which 
has a matrix product form lfT4ll . The number of customers (A^,) decreases approximately 
linearly in time as 

(Nr) ~ (a - J,(p,))t + (No) (6) 

while t < jJ^jLq, starting from a sufficiently large (A^o) at time f = 0. The quantity 
7ii(p,i) is actually the customer current through the right end, i.e. the outflow. The 
system length exhibits a similar behavior 

(L,)^'^-^t.iLo), (7) 

where the density profile is almost flat with the bulk density p,,. In view of the form 
(j5]l, we call the region a < /||(P||) with /3 > /3c "maximal-current-convergent (MC-C) 
phase", and a < Juipn) with /3 < /3c "high-density-convergent (HD-C) phase". 

When a > Juipn), the system does not have a stationary state, and (A^,) and {L,} 
diverge linearly in time. For (A^,), the form (j6|l is valid and we have the asymptotic 
behavior 

{N,)^(a-JM)t (t^oo). (8) 

In view of the form (|5]l, we call the region a > 7||(pii) with /3 > f3c "maximal-current- 
divergent (MC-D) phase", and a > J^ipn) with /3 < /3c "high-density-divergent (HD-D) 
phase". On the other hand, the form (j7]i is not always valid in the divergent phase (see 



eq. (28 i below). The main purpose of this paper is to determine the velocity V^, for 



{Lt) — Vyf as well as the density profile in the divergent phase. 

2.2 Backward sequential update 

We now consider the discrete-time EQP with backward-sequential update (backward 
EQP): first a customer arrives with probability a, and the customer at the right end is 
extracted with probability f3 (if it exists). Then starting from the rightmost particle and 
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going sequentially to the left up to the leftmost particle, we move each particle forward 
with probability p if possible. For example 

101011 1011001 withprobability Q-XyexpxCl -p)x/?x(l -/?), 

111011 ^ mil with probability (1 -a) X (1 -y8) x/?xp x;?, (9) 
-> with probability (1 - a) + a XyS. 

In the first example of transitions (j9]), the customer on the second site can move to the 
rightmost site, thanks to the backward update. On the other hand, he/she cannot move 
in the parallel case, see the first example of Equation Q. 

The current-density relation for the backward-sequential updat^TASEP is 1231 

pp(l -p) 



(10) 



\-pp 

Simulation results imply that the critical line separating the parameter space of the EQP 
into convergent and divergent phases is given by 

' ®^ (Q<l3<Pc\ f ^ (0<^<M 

where p^ = <^ i-VT^ 
(yS. <y6<l), (J3c<p<l), 



with Pc - I - -ijl - p as in the parallel case. The phase diagram is qualitatively similar 
to that of the parallel EQP (|5]l, see Figure|2] In the special case yS = 1 the system always 
has the stationary state P(0) - 1 and f (otherwise) - 0. Therefore we set < yS < 1 in 
the following. 

When a < J^ipJ), we expect the system to converge to a stationary state, and the 
number of customers (A^,) decreases approximately linearly in time as 



{N,) ~{a- JJpS))t + (Nq) 



(12) 



while t < 



(No) 



J (p )^ Starting from a sufficiently large (A^o) at time f = 0. The system 
length also exhibits a similar behavior 



{L,} 



t + (Lo). 



(13) 



When a > J^{pJ), we expect that the system does not has a stationary state and 
(A^,) and {L,) diverge linearly in time. For (A^,), the form ( [T2| is valid and we have the 
asymptotic behavior 

{N,)^{a-JJp^))t (f^<x,). (14) 



On the other hand, for the divergence of the length the form ( 13 1 is not always valid 



(see eq. (34 1 below). 



(11) 



2.3 Limits and special cases 

The discrete-time EQPs have several known models as special cases or limits. The 
following diagram illustrates the relations between the various models: 



Parallel EQP 



p^l 



Rule 184 EQP with 
stochastic boundaries 



As ^ 



As ^ 



Continuous EQP 



Continuous M/M/1 



As ^ 



As 



Backward EQP 



p^l 



Discrete M/M/1 



(15) 



Note that the sitewise and particlewise ordered updates ( 23 1 are identical here. 
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where As is the length of the discrete time step. (A more precise definition of the Hmit 
will be given below.) 



2.3.1 Parallel update with p - I 

The deterministic hopping cases {p - 1) of the discrete-time EQPs correspond to two 
diff'erent processes. The bulk dynamics of the parallel EQP with p = 1 corresponds to 
the rule 184 cellular automaton with stochastic boundaries. It is still an EQP although 
the customer hopping is deterministic Q. The MC-D and MC-C phases vanish in the 
phase diagram (Fig. |2|. This case has an exact "dynamical state" in a matrix product 
form which enables us to derive asymptotic behaviors of the system length, the number 
of customers and the density profile in the limit f — > oo [18 ] : 



Divergent Phase (a > + /?)): 

<L,) = {a-/3+a/3)t+o(t), (N,) = ~^^^""^ f+o(0, P.v/,r 



(-^<1)' 
(x > 1), 
(16) 



Critical Line (a = /3/(l + yS)): 



^Ml+p) 



l+y8 



-erfc 




Pt 



\7r(l+^)^ 



+ o( Vf) , 



(17) 



Convergent Phase (a < /3/(l + /?)): 



a 



, {N,} 



a{l — a) 



where erfc is the complementary error function erfc(x) = e^^' dy. 



, (18) 



2.3.2 Backward-sequential update with p - I 

The backward EQP with p - I is equivalent to the discrete-time M/M/1 queueing 
process which is no longer an EQP as e.g. the state changes to l_^j_l when the 

N N-l 

customer at the server gets service. Thus no empty site appears between the leftmost 
customer and the server, i.e. we have always A^, = L,, if the system starts from the 
empty chain. In the limit f — > oo, the system shows diff'erent behavior, depending on 
the phase: 

• Divergent Phase (a > /T): 

{L,} = {N,} = ia-/3)t + oit), Px,,,^l I [^^JJ' (19) 
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Figure 3: A schematic picture for the density profile in the divergent phase, where x is 
the rescaled position jit. According to the injection probability (rate) a, the rarefaction 
wave is "cut" by the leftmost customer (x - V) and the server (x - 0). 



Critical Line {a — /3): 



Convergent Phase (a < 



Pjt 



erfc 



X 



«(l-y6) V 
Pil-a) 



(20) 



(21) 



2.3.3 Continuous-time update 

Formally, in the continuous-time limit the probabilities a, /3 and p should be replaced 
by aAs + o(As), /3As + o{As) and As + o(As), respectively, and time t is rescaled as 
t/As. Then the continuous-time limits Ai — > of both discrete-time EQPs yield the 
continuous-time EQP studied in |5|. The current-density relation for the continuous- 
time case is simply 

-/contCp) = PP(1 -p), (22) 
and the phase diagram is given as 



C — ■^cont(Pcont) — 



I3(P-P) 
P 

P 

4 



{Q<[i<pl2\ 
(p/2<p<l). 



where pcom = 



1 



(0<l3<p/2), 

(P > p/2). 
(23) 



The continuous -time M/M/1 queueing process is recovered by the continuous -time lim- 
its of the rule 184 case and the discrete-time M/M/1 queue. It is also obtained by the 
p — > oo limit of the continuous-time EQP. 



3 Subphases in the divergent phase 

We consider the TASEP on an infinite chain with the initial densities Pright (at sites 
j < 0) and Pleft (at sites j > 0), where pieft > Pright- When the current J from left to 
right is given by a function of the density p, the rescaled density profile p(x = j/ 1) is 
weU-described by 

f Pright {X < /(Pright)), 
p(x)-j f-\x) (/(Pleft) > X > /(pnghO), (24) 
i Pleft {X > f(pMt)) 

with f(p) = f26l. We will see that, for the EQPs, the density profiles in the 
divergent phase are obtained by cutting this "rarefaction wave" as in Fig. [3] 
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3.1 Parallel case 

From the current-density relation Q for the parallel-update TASEP we have 

p{2p - 1) 



Up) 



.r'(^) 



1 X 



1 



dp Vl -4/?p(l -p) ' 2 l^p(p-x^)' 

We assume that the (rescaled) density profile p „ , (at site xt and time t) has the form 



(25) 



PxtJ 



(x>y,o>x), 
p(x) (y > jc > 0), 



(26) 



where p(x) is given by ( |24[ l with Pright = Pn and pieft = 1. This assumption is supported 
by simulation results. Here V is the velocity of the system length (L,) ^ Vf. 
Under the assumption ([26]l we have 



f(Q'-7™')^f P(X)£/X 



r 

Jo 



(27) 



where both sides are different expressions for the number of customers. Inserting 7°"' 
7ii(Pii) (see Equation (|5]l) into Equation ( 27 1, we find the velocity 



a-J(P|l) 
Pll 



(I), 



lpa-p + l^pa{\-p){\-a) (II), 



(28) 



or 



(HI), 



where 



I 

n 

III 



< v„ < hipi) 



I.e. 



/,(p,)< < /,(!) i.e. MaxJ^^<f|£gj,^^^j<a<p, 



(29) 



i.e. p < < 1. 



Combing this and the form of p,, given in Equation (jSj), we obtain five subphases in the 
divergent phase. In each phase the rescaled density p„ , has a different form (Fig.|4]i: 



HD-D-I: p,„ - 



Plight (y > X > 0), 
{x> y), 

f Pright {Vx> X> 0), 

HD-D-II: p.„,, ^] / '(-x) {y>x>vx\ MC-D-II: p^,,, ^ 

i (.X > y), 

Pright (V'l > X > 0), 
/-'(-«) (v'2>X>Vi). 

1 (y > X > V2), 
(x > y), 



HD-D-III: p,,., 



MC-D-III: p,,,, 



where 

Pright = Pll 



p-2pp+p^ 



f-Hx) (y>x>0), 

ix> V), 

{ f-\x) {v2>x>Q), 
\ 1 (y>x>y2), 

(x > y), 



f-\x)^f-\x). 



(31) 

The HD-D phase is divided into three phases: (I) plateau, (11) plateau-slope and (III) 
plateau-slope-plateau. On the other hand, the MC-D phase is divided into only two 
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phases: (II) slope and (III) slope-plateau. The plateau near the exit does not appear. 
Figures |5] and |6] show simulation results for the velocities and the density profiles, 
respectively, with parameters 
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(32) 
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Figure 4: The subphases of the EQPs; the parallel EQP (top-left), the continuous-time 
EQP (top-right), the backward EQP with ^ < p < I (bottom-left) and the backward 
EQP with < p < ^ (bottom-right). 



3.2 Backward case 



From the current-density relation (lOi for the TASEP with the backward update we 
have 



dJ^ pil-2p + pp^) 

fjx) - — = 

dp (l-ppY 



p p y I + X 



(33) 



As in the parallel-update case, we assume Equation (|26| with p(x) as in ( 24 1, Pright - p. 



and pieft - 1 . From Equation ( 27 1, we find that the velocity for the system length 
{L,) ^ V^t is given by 



2ylp{\ - p)a- pil - a) 
a 



(I) , 

(II) , 

(HI), 



(34) 
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Figure 5: The growth velocities of the system size and the number of customers for 
the parallel (top), backward (middle) and continuous (bottom) EQPs. The simulation 
data were obtained by averaging 10"^ samples. We see that these agree with the lines 
corresponding to (28 1. The parameters are chosen as in Equations (32 1, (37 1 and (44 1. 
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Figure 6: Rescaled density profiles pj, of the parallel (top, t - 10 ), backward (middle, 
t - 10*) and continuous-time (bottom, t = 2500) EQPs. The parameters are chosen as 
in Equations ( 32 1, ( 37 1 and ( 44 1. The simulation data were obtained by averaging 10'* 
samples. 
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where 



I 

II 
III 



< < fJpS) i.e 
/.(pj < < i.e. Max 



P(.p-P) , „ < (y-Pf 

p(l-/3) ^ " - p(l-p)' 



p(i-p)' p ' — 



(35) 



fj l ) > 



I.e. 



< a < 1 , 



When J < p < 1, the case D-III vanishes and the divergent phase is divided into three 
phases (Fig.|4jl. On the other hand, when < p < j, the structure of the subphases 
is qualitatively similar to the parallel case. The density profiles are given by Equation 
(l30ll with 



Plight = p^, vi = fjp^) = 



1-p 



VI = fJl) = 



1-;/ 



(36) 



Figures |5] and |6] show simulation results of the velocities and the density profiles, re- 
spectively, with parameters 



{a,li,p) 



(0.2,0.1,0.36) HD-D-I o (blue), 

(0.4,0.1,0.36) HD-D-I A (red), 

(0.8,0.1,0.36) HD-D-III x (purple), 

(0.3,0.6,0.36) MC-D-II □(orange), 

(0.8,0.6,0.36) MC-D-III -h (green). 



(37) 



3.3 Continuous-time case 

From the current-density relation p2] l for the continuous-time TASEP, we have 

cont 



/cont(p) 



dp 



p{2p - 1). 



(38) 



The velocity ycont of the system length, the subphases and the density profiles can 
be obtained following the same procedure as for the parallel and backward EQPs, or 
simply by taking the continuous-time limits of the results for the two discrete cases: 



yrnnr — 



g-J(Pcunl) _ 

2 ^fpa - p 
a 



'-,a-/3 (I), 



where 



I : < ycont ^ /cont(Pcont) 

II • /cont(Pcont) < ycont ^ /cont(l) 



(II) , 

(III) , 



P(p-P) < < (p-pr 



in : 



/cont(l) ^cont 



Max(^,f)<a<;?, 
i.e. p < a < 1. 



(39) 



(40) 

(41) 
(42) 



The density profiles are given by Equation ( 30 1 with 



1 



Pright = Pcont, Vi = /cont(Pcont) ^ P - 2/3, V2 = /cont(l) ^ P, f \x) ^ fcontix) = 2 2/?' ^^^"^ 

Figures |5] and |6] show simulation results of the velocities and the density profiles, re- 
spectively, with parameters 
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4 On the critical line 

Since (A^,) and (L,) converge to stationary values in the convergent phase, and diverge 
proportional to t in the divergent phase, we expect that they behave as 



{X.} 



0<rA'<l, for X = L,N. 



(45) 
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Figure 7: The exponents yx for the system length X - L and the number of particles 
on the critical lines of the parallel (p = 0.84, Pc = 0.6, left) and backward (p = 
0.36, /?c = 0.2 and p = 0.19, = 0.1, right) EQPs. The markers □ (orange) and x 
(purple) correspond to ji and j^, respectively, for the parallel case with p = 0.84 and 
the backward case with p - 0.36, and o (blue) and + (red) correspond to ji and y^i, 
respectively, for the backward case with p = 0.19. 



Under this assumption we have 

InX, - \nX,/t 



Jx (t^oc), (46) 



Inb 

for X = L,N. After verifying that the growth behavior is indeed well-described by 



power-laws of the form (45 i, we estimate the exponents yx by applying (46 1 to simu- 
lated samples with b = 10 and f < 5 x 10^. The number of samples for this estimation 
for each parameter set is basically 10"^, but 10^ or 5 x lO'' samples were used for the 
backward EQP in the region 0.28 < P < 0.8 and (3 > 0.8, respectively, because there 
fluctuations of L and are very large. The results shown in Fig. [7] are consistent with 
the expectation - y^ everywhere on the critical line. This is supported by the ob- 
servation that the total density ptot - {Ni}/{L,} reaches quickly an almost stationary 
value which implies that yi = y^. More detailed results will be present in a future 
publication. 

The critical lines of the EQPs consist of two parts: a curved and a straight line 
(Fig.|2]i. On the curved part, the simulation results indicate 

<M),<L,)-0(Vf) . (47) 
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The behavior of <L,) and (A^,) on the straight part of the critical line is not so clear 
although diffusive behavior can be excluded. As Fig. |7] indicates, the exponents are 
smaller than on the curved part, i.e. Jl = 7n < ^ /2. For the parallel case, y - 1/4 (with 
large corrections near /3c) can not be excluded, but for the backward case, the exponents 
seems to depend on the value of p. For example, the exponents for p = 0.19 seem to 
be bigger than those for p - 0.36, see the right graph of Fig. |7] Our simulation results 
are not sufficient to determine conclusively the dependency of y on the parameters, e.g. 
how it varies with /3 near /3c. 



5 Conclusion 

We have continued our studies of the exclusive queueing process (EQP) which extends 
the classical M/M/1 queueing process by incorporating the excluded volume effect. 
We have compared the behavior of the model with different update schemes (parallel, 
backward-sequential, continuous time). The phase diagrams are qualitatively similar, 
except for certain limiting cases. 

The phase diagram of the EQP turns out to be rather rich. Here we have shown 
that the divergent phase is subdivided into up to 5 different subphases according to the 
parameter dependence of the current and the density profiles. The MC-D phase has 
two different subphases (the slope and plateau-slope phases), and the HD-D phase has 
three different subphases (the plateau, plateau-slope and plateau-slope-plateau phases). 

In the divergent phase we have conjectured the analytic form of the density profiles 
which show good agreement with simulation results. The shapes of the rescaled profiles 
can be understood in terms of a rarefaction wave that is "cut" at both ends. 

On the critical line separating the divergent from the convergent phase the length of 
the system grows sublinearly. Based on simulation results we find diffusive behavior 
on the curved part of the critical line (i.e. /3 < /3c) for all updates. In the special 
case p = I for the two discrete EQPs, the density profiles can be written in terms of 



the complementary error function as Equations ( 17i and (20 1. Identifying the density 
profile on the curved part for the EQPs with general values of p is one of problems that 
need to be clarified in the future. 

The behavior on the straight part /? > /3c of the critical line is subdiffusive (y < 
1/2). However we could not clearly determine whether the exponent y depends on the 
parameters, and more simulation data with sufficient accuracy are needed to determine 
the behavior on the straight part. 
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